rm(list=c(ls()))

library(rio)
library(ggplot2)
library(stargazer)
library(foreign)
library(margins)
library(interplot)

interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

data<-import("data_Bertelsmann_2020-03-short.xlsx")

#coding--------------
#female
data$female<-ifelse(data$gender=="female", 1,0)

#education
data$education<- car::recode(data$education, "2=1; 3=2; 4=3")

#unemployment
data$unemployed<-ifelse(data$employment==2, 1,0)

#italian clusre
data$ITclosure <- ifelse(data$day>8, 1,0)
table(data$pmpar)
#PM Partisans
data$pmpar=0
data$pmpar[data$party_fr=="RenM" | data$party_de=="CDU-CSU" | data$party_pl=="PiS" |data$party_es=="PSOE"]=1

#Government partisans
data$govpar=0
data$govpar[data$party_fr=="RenM" | data$party_de=="CDU-CSU" | data$party_de=="SPD" | data$party_pl=="PiS" |data$party_es=="PSOE" | data$party_es=="Podemos" ]=1

#Populist right
data$rpppar=0
data$rpppar[data$party_fr=="RN" | data$party_de=="AfD" | data$party_pl=="PiS" |data$party_es=="Vox"]=1

#left-right
data$lr<-car::recode(data$leftright, "2=1; 3=2; 4=3; 5=4; 6=4")

#placebo treatments
data$placebo_closure_march5 <-ifelse(data$day==5, 0,1)
data$placebo_closure_march6 <-ifelse(data$day<=6, 0,1)
data$placebo_closure_march7 <-ifelse(data$day<=7, 0,1)


#select data-------------------------
#Full models
cols<- c("govpar", "pmpar", "rpppar", "lr", "ITclosure", "female", "age", "education", "class", "city", "unemployed", "country_code", "weight")
datasub <- data[cols]


#Table 1------------------
summary(lm(datasub$age ~ datasub$ITclosure, weights=datasub$weights))
summary(lm(datasub$female ~ datasub$ITclosure, weights=datasub$weights))
summary(lm(datasub$city ~ datasub$ITclosure, weights=datasub$weights))
summary(lm(datasub$education ~ datasub$ITclosure, weights=datasub$weights))
summary(lm(datasub$class ~ datasub$ITclosure, weights=datasub$weights))
summary(lm(datasub$unemployed ~ datasub$ITclosure, weights=datasub$weights))

#Figure 3 in the main text-----------------------

#outcome
options <- expand.grid(colnames(datasub)[1:4],colnames(datasub)[5]) # 1-4 are different DVs, 5 = IV
data_catherine <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, weights=weight)
  data_catherine[i,1] <- coef(model)[2]
  data_catherine[i,2] <- summary(model)$coefficients[2,2]
}

#Models withonly country fixed effects
options <- expand.grid(colnames(datasub)[1:4],colnames(datasub)[5]) # 5&6 are DVs, 1-4 are different IVs
data_basic <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) +  as.factor(country_code), data=datasub, weights=weight)
  data_basic[i,1] <- coef(model)[2]
  data_basic[i,2] <- summary(model)$coefficients[2,2]
}

#build plot
forplot<- rbind(data_catherine, data_basic)
forplot$dv<-rep(c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"),2)
forplot$dv <- factor(forplot$dv,levels = c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"))
forplot$controls<-c(rep("Yes",4), rep("No", 4))

main<-ggplot(forplot, aes(x = controls)) + facet_grid(.~dv, scales = "free_y", space="free")+  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +  geom_pointrange(aes(x = controls, colour=controls, y = X1, ymin = X1 - X2*interval2, ymax = X1 + X2*interval2), lwd = 1/2, position = position_dodge(width = 1/2), shape = 21) + theme_bw() + xlab("Control variables included (Yes/No)") + ylab("Unstandardized OLS regression coefficient")+ geom_linerange(aes(x = controls, colour=controls, ymin = X1 - X2*interval1,  ymax = X1 + X2*interval1), lwd = 1, position = position_dodge(width = 1/2))+ scale_color_manual(values=c('Grey','Black'), guide = FALSE)
ggsave(main, file="Main_results.pdf")

#Figure 4: restrict sample to day when national lockdown was enforced in country ------------
#code country
data$ITclosure_stop_country_lock<-NA
data$ITclosure_stop_country_lock[data$day<9 & data$country_code=="PL" | data$day<9 & data$country_code=="ES" | data$day<9 & data$country_code=="FR" | data$day<9 & data$country_code=="DE"]=0
data$ITclosure_stop_country_lock[data$day>8 & data$day<13 & data$country_code=="PL" | data$day>8 & data$day<14 & data$country_code=="ES" | data$day>8 & data$day<17 & data$country_code=="FR"| data$day>8 & data$day<23 & data$country_code=="DE"]=1

cols<- c("govpar", "pmpar", "rpppar", "lr", "ITclosure_stop_country_lock", "female", "age", "education", "class", "unemployed", "country_code", "weight")

datarobust <- data[cols]

options <- expand.grid(colnames(datarobust)[1:4],colnames(datarobust)[5]) # 1-4 are different DVs, 5 = IV
data_catherine <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) + female + age + education + class + unemployed +  as.factor(country_code), data=datarobust, weights=weight)
  data_catherine[i,1] <- coef(model)[2]
  data_catherine[i,2] <- summary(model)$coefficients[2,2]
}

#Models withonly country fixed effects
options <- expand.grid(colnames(datarobust)[1:4],colnames(datarobust)[5]) # 5&6 are DVs, 1-4 are different IVs
data_basic <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) +  as.factor(country_code), data=datarobust, weights=weight)
  data_basic[i,1] <- coef(model)[2]
  data_basic[i,2] <- summary(model)$coefficients[2,2]
}


forplot<- rbind(data_catherine, data_basic)
forplot$dv<-rep(c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"),2)
forplot$dv <- factor(forplot$dv,levels = c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"))
forplot$controls<-c(rep("Yes",4), rep("No", 4))

main<-ggplot(forplot, aes(x = controls)) + facet_grid(.~dv, scales = "free_y", space="free")+  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +  geom_pointrange(aes(x = controls, colour=controls, y = X1, ymin = X1 - X2*interval2, ymax = X1 + X2*interval2), lwd = 1/2, position = position_dodge(width = 1/2), shape = 21) + theme_bw() + xlab("Control variables included (Yes/No)") + ylab("Unstandardized OLS regression coefficient")+ geom_linerange(aes(x = controls, colour=controls, ymin = X1 - X2*interval1,  ymax = X1 + X2*interval1), lwd = 1, position = position_dodge(width = 1/2))+ scale_color_manual(values=c('Grey','Black'), guide = FALSE)
ggsave(main, file="Robustness_main_constrained_lockdown.pdf")






#Figure 5: Placebo test ------------
#Full models
cols<- c("govpar", "pmpar", "rpppar", "placebo_closure_march5", "placebo_closure_march6", "placebo_closure_march7", "gender", "age", "city", "education", "class", "unemployed", "country_code", "weight")
data_placebo <- data[cols]

#outcome
options <- expand.grid(colnames(data_placebo)[1:3],colnames(data_placebo)[4:6]) # 5&6 are DVs, 1-4 are different IVs
data_catherine <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) + gender + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=data_placebo, weights=weight)
  data_catherine[i,1] <- coef(model)[2]
  data_catherine[i,2] <- summary(model)$coefficients[2,2]
}

#Models withonly country fixed effects
options <- expand.grid(colnames(data_placebo)[1:3],colnames(data_placebo)[4:6]) # 5&6 are DVs, 1-4 are different IVs
data_basic <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) +  as.factor(country_code), data=data_placebo, weights=weight)
  data_basic[i,1] <- coef(model)[2]
  data_basic[i,2] <- summary(model)$coefficients[2,2]
}

forplot<- rbind(data_catherine, data_basic)
forplot$dv<-rep(c("Government Party\n support", "PM Party\n support", "Far Right Party\n support"),6)
forplot$dv <- factor(forplot$dv,levels = c("Government Party\n support", "PM Party\n support", "Far Right Party\n support"))


forplot$controls<-c(rep("Yes",9), rep("No", 9))
forplot$data<-c(rep("March 5",3), rep("March 6", 3), rep("March 7", 3))
forplot$data <- factor(forplot$data,levels = c("March 5", "March 6", "March 7"))

forplot$grouping_colour <-ifelse(forplot$controls=="yes", 1,0)
placebo<-ggplot(forplot, aes(x = controls)) + facet_grid(data~dv, scales = "free_y", space="free")+  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +  geom_pointrange(aes(x = controls, colour=controls, y = X1, ymin = X1 - X2*interval2, ymax = X1 + X2*interval2), lwd = 1/2, position = position_dodge(width = 1/2), shape = 21) + theme_bw() + xlab("Control variables included (Yes/No)") + ylab("Unstandardized OLS regression coefficient")+ geom_linerange(aes(x = controls, colour=controls, ymin = X1 - X2*interval1,  ymax = X1 + X2*interval1), lwd = 1, position = position_dodge(width = 1/2))+ scale_color_manual(values=c('Grey','Black'), guide = FALSE)
ggsave(placebo, file="placebo_results.pdf")








#Figure 6: Conditional upon education-----------------
dataint<-datasub
dataint$education<-as.factor(dataint$education)
interaction<-list()
interaction[[1]] <- lm(govpar ~ ITclosure*education + female + age + city +  + as.factor(class) + unemployed +  as.factor(country_code), data=dataint, weights=weight)
interaction[[2]] <- lm(pmpar ~ ITclosure*education + female + age + city +  + as.factor(class) + unemployed +  as.factor(country_code), data=dataint, weights=weight)
interaction[[3]] <- lm(rpppar ~ ITclosure*education + female + age + city +  + as.factor(class) + unemployed +  as.factor(country_code), data=dataint, weights=weight)
interaction[[4]] <- lm(lr ~ ITclosure*education + female + age + city +  + as.factor(class) + unemployed +  as.factor(country_code), data=dataint, weights=weight)


a <- summary(margins(interaction[[1]],dataint,variables ="ITclosure" ,at=list(education=c("1","2","3"))),level=.95)
b <- summary(margins(interaction[[2]],dataint,variables ="ITclosure" ,at=list(education=c("1","2","3"))),level=.95)
c <- summary(margins(interaction[[3]],dataint,variables ="ITclosure" ,at=list(education=c("1","2","3"))),level=.95)
d <- summary(margins(interaction[[3]],dataint,variables ="ITclosure" ,at=list(education=c("1","2","3"))),level=.95)

out <- rbind(a,b,c,d)
out$education<-as.factor(out$education)
levels(out$education) <- c("Some high school", "High school", "University")
out$split <- c(1,1,1,2,2,2,3,3,3,4,4,4)
out$split<-as.factor(out$split)
levels(out$split) <- c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology")

edu_plot<-ggplot(out,aes(x=education,y=AME))+geom_pointrange(aes(ymin=lower,ymax=upper))+facet_wrap(~split, ncol=2)+theme_bw()+ xlab("Educational levels") + ylab("Average Marginal Effect of Italian lockdown on dependent variables") + geom_hline(yintercept=0, linetype="dashed", color = "red")
ggsave(edu_plot, file="conditional_education.pdf")



#Appendix Table A.1 and A.2 - Results belonging to Figure 3 in appendix---------------
main_results<-list()
main_results[[1]] <- lm(govpar ~ ITclosure + as.factor(country_code), data=datasub, weights=weight)
main_results[[2]] <- lm(govpar ~ ITclosure + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, weights=weight)
main_results[[3]] <- lm(pmpar ~ ITclosure + as.factor(country_code), data=datasub, weights=weight)
main_results[[4]] <- lm(pmpar ~ ITclosure + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, weights=weight)
main_results[[5]] <- lm(rpppar ~ ITclosure + as.factor(country_code), data=datasub, weights=weight)
main_results[[6]] <- lm(rpppar ~ ITclosure + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, weights=weight)
main_results[[7]] <- lm(lr ~ ITclosure + as.factor(country_code), data=datasub, weights=weight)
main_results[[8]] <- lm(lr ~ ITclosure + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, weights=weight)

stargazer(main_results[[1]], main_results[[3]], main_results[[5]], main_results[[7]], column.separate = c(1,1,1,1), dep.var.caption  = "Models without control variables", dep.var.labels = c("Government Party", "PM Party", "Far Right Party", "Left-right"), covariate.labels = c("Lockdown", "Country: Spain", "Country: France", "Country: Poland"), title="Results belonging to Figure 3 in the main text: models without control variables", no.space = TRUE, align=TRUE, model.numbers = FALSE, omit.stat=c("LL","ser","f", "adj.rsq"), star.char = c("+", "*"), star.cutoffs = c(.1, 0.05), notes = c("+ p<.1; * p<0.05)"), notes.append = F, out="Fig2_base.tex", label="tab:fig2_base", digits=3)

stargazer(main_results[[2]], main_results[[4]], main_results[[6]], main_results[[8]], column.separate = c(1,1,1,1), dep.var.caption  = "Models with control variables", dep.var.labels = c("Government Party", "PM Party", "Far Right Party", "Left-right"), covariate.labels = c("Lockdown", "Female", "Age", "Urban", "Education: high school", "Education: university", "Class: lower middle", "Class: upper middle", "Class: upper class", "Unemployed", "Country: Spain", "Country: France", "Country: Poland"), title="Results belonging to Figure 3 in the main text: models without control variables", no.space = TRUE, align=TRUE, model.numbers = FALSE, omit.stat=c("LL","ser","f", "adj.rsq"), star.char = c("+", "*"), star.cutoffs = c(.1, 0.05), notes = c("+ p<.1; * p<0.05)"), notes.append = F, out="Fig2_full.tex", label="tab:fig2_full", digits=3)


#Appendix Figure A.1 - Results belonging to restircting it to 1 week after the lockdown-------------------------------------
data$ITclosure_1week<-NA
data$ITclosure_1week[data$day<9]=0
data$ITclosure_1week[data$day>8 & data$day<16]=1

cols<- c("govpar", "pmpar", "rpppar", "lr", "ITclosure_1week", "female", "age", "education", "class", "unemployed", "country_code", "weight")
datarobust <- data[cols]

options <- expand.grid(colnames(datarobust)[1:4],colnames(datarobust)[5]) # 1-4 are different DVs, 5 = IV
data_catherine <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) + female + age + education + class + unemployed +  as.factor(country_code), data=datarobust, weights=weight)
  data_catherine[i,1] <- coef(model)[2]
  data_catherine[i,2] <- summary(model)$coefficients[2,2]
}

#Models withonly country fixed effects
options <- expand.grid(colnames(datarobust)[1:4],colnames(datarobust)[5]) # 5&6 are DVs, 1-4 are different IVs
data_basic <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) +  as.factor(country_code), data=datarobust, weights=weight)
  data_basic[i,1] <- coef(model)[2]
  data_basic[i,2] <- summary(model)$coefficients[2,2]
}

forplot<- rbind(data_catherine, data_basic)
forplot$dv<-rep(c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"),2)
forplot$dv <- factor(forplot$dv,levels = c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"))
forplot$controls<-c(rep("Yes",4), rep("No", 4))

main<-ggplot(forplot, aes(x = controls)) + facet_grid(.~dv, scales = "free_y", space="free")+  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +  geom_pointrange(aes(x = controls, colour=controls, y = X1, ymin = X1 - X2*interval2, ymax = X1 + X2*interval2), lwd = 1/2, position = position_dodge(width = 1/2), shape = 21) + theme_bw() + xlab("Control variables included (Yes/No)") + ylab("Unstandardized OLS regression coefficient")+ geom_linerange(aes(x = controls, colour=controls, ymin = X1 - X2*interval1,  ymax = X1 + X2*interval1), lwd = 1, position = position_dodge(width = 1/2))+ scale_color_manual(values=c('Grey','Black'), guide = FALSE)
ggsave(main, file="Robustness_main_1week.pdf")







#Appendix Figure A.2 - control for deaths or cases at the time -----------------------------
data$deaths_at_survey<-NA
data$deaths_at_survey[data$day==5 & data$country_code=="PL"]=0
data$deaths_at_survey[data$day==6 & data$country_code=="PL"]=0
data$deaths_at_survey[data$day==7 & data$country_code=="PL"]=0
data$deaths_at_survey[data$day==8 & data$country_code=="PL"]=0
data$deaths_at_survey[data$day==9 & data$country_code=="PL"]=0
data$deaths_at_survey[data$day==10 & data$country_code=="PL"]=0
data$deaths_at_survey[data$day==11 & data$country_code=="PL"]=0
data$deaths_at_survey[data$day==12 & data$country_code=="PL"]=1
data$deaths_at_survey[data$day==13 & data$country_code=="PL"]=1
data$deaths_at_survey[data$day==14 & data$country_code=="PL"]=1
data$deaths_at_survey[data$day==15 & data$country_code=="PL"]=3
data$deaths_at_survey[data$day==16 & data$country_code=="PL"]=3
data$deaths_at_survey[data$day==17 & data$country_code=="PL"]=5
data$deaths_at_survey[data$day==18 & data$country_code=="PL"]=5
data$deaths_at_survey[data$day==19 & data$country_code=="PL"]=5
data$deaths_at_survey[data$day==20 & data$country_code=="PL"]=5
data$deaths_at_survey[data$day==21 & data$country_code=="PL"]=5
data$deaths_at_survey[data$day==22 & data$country_code=="PL"]=5
data$deaths_at_survey[data$day==23 & data$country_code=="PL"]=7
data$deaths_at_survey[data$day==24 & data$country_code=="PL"]=8
data$deaths_at_survey[data$day==25 & data$country_code=="PL"]=10

#Germany
data$deaths_at_survey[data$day==5 & data$country_code=="DE"]=0
data$deaths_at_survey[data$day==6 & data$country_code=="DE"]=0
data$deaths_at_survey[data$day==7 & data$country_code=="DE"]=0
data$deaths_at_survey[data$day==8 & data$country_code=="DE"]=0
data$deaths_at_survey[data$day==9 & data$country_code=="DE"]=2
data$deaths_at_survey[data$day==10 & data$country_code=="DE"]=2
data$deaths_at_survey[data$day==11 & data$country_code=="DE"]=3
data$deaths_at_survey[data$day==12 & data$country_code=="DE"]=3
data$deaths_at_survey[data$day==13 & data$country_code=="DE"]=3
data$deaths_at_survey[data$day==14 & data$country_code=="DE"]=8
data$deaths_at_survey[data$day==15 & data$country_code=="DE"]=12
data$deaths_at_survey[data$day==16 & data$country_code=="DE"]=12
data$deaths_at_survey[data$day==17 & data$country_code=="DE"]=12
data$deaths_at_survey[data$day==18 & data$country_code=="DE"]=12
data$deaths_at_survey[data$day==19 & data$country_code=="DE"]=20
data$deaths_at_survey[data$day==20 & data$country_code=="DE"]=45
data$deaths_at_survey[data$day==21 & data$country_code=="DE"]=67
data$deaths_at_survey[data$day==22 & data$country_code=="DE"]=67
data$deaths_at_survey[data$day==23 & data$country_code=="DE"]=94
data$deaths_at_survey[data$day==24 & data$country_code=="DE"]=126
data$deaths_at_survey[data$day==25 & data$country_code=="DE"]=149

#Spain
data$deaths_at_survey[data$day==5 & data$country_code=="ES"]=3
data$deaths_at_survey[data$day==6 & data$country_code=="ES"]=5
data$deaths_at_survey[data$day==7 & data$country_code=="ES"]=5
data$deaths_at_survey[data$day==8 & data$country_code=="ES"]=10
data$deaths_at_survey[data$day==9 & data$country_code=="ES"]=28
data$deaths_at_survey[data$day==10 & data$country_code=="ES"]=36
data$deaths_at_survey[data$day==11 & data$country_code=="ES"]=47
data$deaths_at_survey[data$day==12 & data$country_code=="ES"]=84
data$deaths_at_survey[data$day==13 & data$country_code=="ES"]=121
data$deaths_at_survey[data$day==14 & data$country_code=="ES"]=137
data$deaths_at_survey[data$day==15 & data$country_code=="ES"]=289
data$deaths_at_survey[data$day==16 & data$country_code=="ES"]=310
data$deaths_at_survey[data$day==17 & data$country_code=="ES"]=492
data$deaths_at_survey[data$day==18 & data$country_code=="ES"]=599
data$deaths_at_survey[data$day==19 & data$country_code=="ES"]=768
data$deaths_at_survey[data$day==20 & data$country_code=="ES"]=1003
data$deaths_at_survey[data$day==21 & data$country_code=="ES"]=1327
data$deaths_at_survey[data$day==22 & data$country_code=="ES"]=1327
data$deaths_at_survey[data$day==23 & data$country_code=="ES"]=1721
data$deaths_at_survey[data$day==24 & data$country_code=="ES"]=2183
data$deaths_at_survey[data$day==25 & data$country_code=="ES"]=2697

#France
data$deaths_at_survey[data$day==5 & data$country_code=="FR"]=7
data$deaths_at_survey[data$day==6 & data$country_code=="FR"]=10
data$deaths_at_survey[data$day==7 & data$country_code=="FR"]=11
data$deaths_at_survey[data$day==8 & data$country_code=="FR"]=20
data$deaths_at_survey[data$day==9 & data$country_code=="FR"]=31
data$deaths_at_survey[data$day==10 & data$country_code=="FR"]=34
data$deaths_at_survey[data$day==11 & data$country_code=="FR"]=49
data$deaths_at_survey[data$day==12 & data$country_code=="FR"]=49
data$deaths_at_survey[data$day==13 & data$country_code=="FR"]=80
data$deaths_at_survey[data$day==14 & data$country_code=="FR"]=92
data$deaths_at_survey[data$day==15 & data$country_code=="FR"]=128
data$deaths_at_survey[data$day==16 & data$country_code=="FR"]=149
data$deaths_at_survey[data$day==17 & data$country_code=="FR"]=176
data$deaths_at_survey[data$day==18 & data$country_code=="FR"]=245
data$deaths_at_survey[data$day==19 & data$country_code=="FR"]=373
data$deaths_at_survey[data$day==20 & data$country_code=="FR"]=451
data$deaths_at_survey[data$day==21 & data$country_code=="FR"]=563
data$deaths_at_survey[data$day==22 & data$country_code=="FR"]=563
data$deaths_at_survey[data$day==23 & data$country_code=="FR"]=675
data$deaths_at_survey[data$day==24 & data$country_code=="FR"]=861
data$deaths_at_survey[data$day==25 & data$country_code=="FR"]=1101

#robustness test 2: control for case at the time
data$cases_at_survey<-NA
data$cases_at_survey[data$day==5 & data$country_code=="PL"]=1
data$cases_at_survey[data$day==6 & data$country_code=="PL"]=5
data$cases_at_survey[data$day==7 & data$country_code=="PL"]=6
data$cases_at_survey[data$day==8 & data$country_code=="PL"]=11
data$cases_at_survey[data$day==9 & data$country_code=="PL"]=16
data$cases_at_survey[data$day==10 & data$country_code=="PL"]=22
data$cases_at_survey[data$day==11 & data$country_code=="PL"]=27
data$cases_at_survey[data$day==12 & data$country_code=="PL"]=51
data$cases_at_survey[data$day==13 & data$country_code=="PL"]=64
data$cases_at_survey[data$day==14 & data$country_code=="PL"]=64
data$cases_at_survey[data$day==15 & data$country_code=="PL"]=111
data$cases_at_survey[data$day==16 & data$country_code=="PL"]=150
data$cases_at_survey[data$day==17 & data$country_code=="PL"]=246
data$cases_at_survey[data$day==18 & data$country_code=="PL"]=287
data$cases_at_survey[data$day==19 & data$country_code=="PL"]=325
data$cases_at_survey[data$day==20 & data$country_code=="PL"]=425
data$cases_at_survey[data$day==21 & data$country_code=="PL"]=536
data$cases_at_survey[data$day==22 & data$country_code=="PL"]=634
data$cases_at_survey[data$day==23 & data$country_code=="PL"]=749
data$cases_at_survey[data$day==24 & data$country_code=="PL"]=901
data$cases_at_survey[data$day==25 & data$country_code=="PL"]=1051

#Germany
data$cases_at_survey[data$day==5 & data$country_code=="DE"]=400
data$cases_at_survey[data$day==6 & data$country_code=="DE"]=639
data$cases_at_survey[data$day==7 & data$country_code=="DE"]=795
data$cases_at_survey[data$day==8 & data$country_code=="DE"]=902
data$cases_at_survey[data$day==9 & data$country_code=="DE"]=1139
data$cases_at_survey[data$day==10 & data$country_code=="DE"]=1296
data$cases_at_survey[data$day==11 & data$country_code=="DE"]=1567
data$cases_at_survey[data$day==12 & data$country_code=="DE"]=1577
data$cases_at_survey[data$day==13 & data$country_code=="DE"]=2260
data$cases_at_survey[data$day==14 & data$country_code=="DE"]=3795
data$cases_at_survey[data$day==15 & data$country_code=="DE"]=4838
data$cases_at_survey[data$day==16 & data$country_code=="DE"]=6012
data$cases_at_survey[data$day==17 & data$country_code=="DE"]=7156
data$cases_at_survey[data$day==18 & data$country_code=="DE"]=8198
data$cases_at_survey[data$day==19 & data$country_code=="DE"]=10999
data$cases_at_survey[data$day==20 & data$country_code=="DE"]=18323
data$cases_at_survey[data$day==21 & data$country_code=="DE"]=21463
data$cases_at_survey[data$day==22 & data$country_code=="DE"]=21463
data$cases_at_survey[data$day==23 & data$country_code=="DE"]=24774
data$cases_at_survey[data$day==24 & data$country_code=="DE"]=29212
data$cases_at_survey[data$day==25 & data$country_code=="DE"]=31554

#Spain
data$cases_at_survey[data$day==5 & data$country_code=="ES"]=257
data$cases_at_survey[data$day==6 & data$country_code=="ES"]=374
data$cases_at_survey[data$day==7 & data$country_code=="ES"]=430
data$cases_at_survey[data$day==8 & data$country_code=="ES"]=589
data$cases_at_survey[data$day==9 & data$country_code=="ES"]=1024
data$cases_at_survey[data$day==10 & data$country_code=="ES"]=1639
data$cases_at_survey[data$day==11 & data$country_code=="ES"]=2140
data$cases_at_survey[data$day==12 & data$country_code=="ES"]=2965
data$cases_at_survey[data$day==13 & data$country_code=="ES"]=4231
data$cases_at_survey[data$day==14 & data$country_code=="ES"]=5753
data$cases_at_survey[data$day==15 & data$country_code=="ES"]=7753
data$cases_at_survey[data$day==16 & data$country_code=="ES"]=9191
data$cases_at_survey[data$day==17 & data$country_code=="ES"]=11178
data$cases_at_survey[data$day==18 & data$country_code=="ES"]=13716
data$cases_at_survey[data$day==19 & data$country_code=="ES"]=17147
data$cases_at_survey[data$day==20 & data$country_code=="ES"]=19980
data$cases_at_survey[data$day==21 & data$country_code=="ES"]=24926
data$cases_at_survey[data$day==22 & data$country_code=="ES"]=24926
data$cases_at_survey[data$day==23 & data$country_code=="ES"]=28572
data$cases_at_survey[data$day==24 & data$country_code=="ES"]=33089
data$cases_at_survey[data$day==25 & data$country_code=="ES"]=39673

#France
data$cases_at_survey[data$day==5 & data$country_code=="FR"]=420
data$cases_at_survey[data$day==6 & data$country_code=="FR"]=613
data$cases_at_survey[data$day==7 & data$country_code=="FR"]=706
data$cases_at_survey[data$day==8 & data$country_code=="FR"]=1116
data$cases_at_survey[data$day==9 & data$country_code=="FR"]=1402
data$cases_at_survey[data$day==10 & data$country_code=="FR"]=1774
data$cases_at_survey[data$day==11 & data$country_code=="FR"]=2269
data$cases_at_survey[data$day==12 & data$country_code=="FR"]=2281
data$cases_at_survey[data$day==13 & data$country_code=="FR"]=3640
data$cases_at_survey[data$day==14 & data$country_code=="FR"]=4469
data$cases_at_survey[data$day==15 & data$country_code=="FR"]=5380
data$cases_at_survey[data$day==16 & data$country_code=="FR"]=6573
data$cases_at_survey[data$day==17 & data$country_code=="FR"]=7652
data$cases_at_survey[data$day==18 & data$country_code=="FR"]=9043
data$cases_at_survey[data$day==19 & data$country_code=="FR"]=10877
data$cases_at_survey[data$day==20 & data$country_code=="FR"]=12475
data$cases_at_survey[data$day==21 & data$country_code=="FR"]=14296
data$cases_at_survey[data$day==22 & data$country_code=="FR"]=14296
data$cases_at_survey[data$day==23 & data$country_code=="FR"]=15821
data$cases_at_survey[data$day==24 & data$country_code=="FR"]=19615
data$cases_at_survey[data$day==25 & data$country_code=="FR"]=22025

#select data
cols<- c("govpar", "pmpar", "rpppar", "lr", "ITclosure", "female", "age", "education", "class", "unemployed", "country_code", "weight", "cases_at_survey", "deaths_at_survey")
datacases <- data[cols]

#control for deaths at time of the survey
options <- expand.grid(colnames(datasub)[1:4],colnames(datasub)[5]) # 1-4 are different DVs, 5 = IV
data_catherine <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) + female + age + education + class + unemployed +  as.factor(country_code) + deaths_at_survey, data=datacases, weights=weight)
  data_catherine[i,1] <- coef(model)[2]
  data_catherine[i,2] <- summary(model)$coefficients[2,2]
}

#Models withonly country fixed effects
options <- expand.grid(colnames(datasub)[1:4],colnames(datasub)[5]) # 5&6 are DVs, 1-4 are different IVs
data_basic <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) + female + age + education + class + unemployed +  as.factor(country_code) + cases_at_survey, data=datacases, weights=weight)
  data_basic[i,1] <- coef(model)[2]
  data_basic[i,2] <- summary(model)$coefficients[2,2]
}

#build plot
forplot<- rbind(data_catherine, data_basic)
forplot$dv<-rep(c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"),2)
forplot$dv <- factor(forplot$dv,levels = c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"))
forplot$controls<-c(rep("Deaths",4), rep("Cases", 4))


main<-ggplot(forplot, aes(x = controls)) + facet_grid(.~dv, scales = "free_y", space="free")+  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +  geom_pointrange(aes(x = controls, colour=controls, y = X1, ymin = X1 - X2*interval2, ymax = X1 + X2*interval2), lwd = 1/2, position = position_dodge(width = 1/2), shape = 21) + theme_bw() + xlab("Control variables included (Yes/No)") + ylab("Unstandardized OLS regression coefficient")+ geom_linerange(aes(x = controls, colour=controls, ymin = X1 - X2*interval1,  ymax = X1 + X2*interval1), lwd = 1, position = position_dodge(width = 1/2))+ scale_color_manual(values=c('Grey','Black'), guide = FALSE)
ggsave(main, file="Robustness_control_cases_deaths.pdf")

#Appendix Table A.3 - logistic regrestic ------------------------
log_results<-list()
log_results[[1]] <- glm(govpar ~ ITclosure + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, family = binomial, weights=weight)
log_results[[2]] <- glm(pmpar ~ ITclosure + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, family = binomial, weights=weight)
log_results[[3]] <- glm(rpppar ~ ITclosure + female + age + city + as.factor(education) + as.factor(class) + unemployed +  as.factor(country_code), data=datasub, family = binomial, weights=weight)

stargazer(log_results[[1]], log_results[[2]], log_results[[3]], column.separate = c(1,1,1), dep.var.caption  = "Models with control variables", dep.var.labels = c("Government Party", "PM Party", "Far Right Party"), covariate.labels = c("Lockdown", "Female", "Age", "Urban", "Education: high school", "Education: university", "Class: lower middle", "Class: upper middle", "Class: upper class", "Unemployed", "Country: Spain", "Country: France", "Country: Poland"), title="Results belonging to Figure 3 in the main text: logistic regression results", no.space = TRUE, align=TRUE, model.numbers = FALSE, omit.stat=c("LL","ser","f", "adj.rsq"), star.char = c("+", "*"), star.cutoffs = c(.1, 0.05), notes = c("+ p<.1; * p<0.05)"), notes.append = F, out="Fig2_log.tex", label="tab:fig2_log", digits=3)


#Appendix Figure A.3 - without survey weights------------------
options <- expand.grid(colnames(datasub)[1:4],colnames(datasub)[5]) # 5&6 are DVs, 1-4 are different IVs
data_catherine <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) + female + age + education + class + unemployed +  as.factor(country_code), data=datasub)
  data_catherine[i,1] <- coef(model)[2]
  data_catherine[i,2] <- summary(model)$coefficients[2,2]
}

#Models withonly country fixed effects
options <- expand.grid(colnames(datasub)[1:4],colnames(datasub)[5]) # 5&6 are DVs, 1-4 are different IVs
data_basic <- data.frame(matrix(NA,dim(options)[1],2),options)

#loop
for (i in 1:dim(options)[1]){
  outcome <- as.character(options[i,1])
  exposure <- as.character(options[i,2])
  model <- lm(get(outcome) ~ get(exposure) +  as.factor(country_code), data=datasub)
  data_basic[i,1] <- coef(model)[2]
  data_basic[i,2] <- summary(model)$coefficients[2,2]
}

#build plot
forplot<- rbind(data_catherine, data_basic)
forplot$dv<-rep(c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"),2)
forplot$dv <- factor(forplot$dv,levels = c("Government Party\n support", "PM Party\n support", "Far Right Party\n support", "Left-right\n ideology"))

forplot$controls<-c(rep("Yes",4), rep("No", 4))

main<-ggplot(forplot, aes(x = controls)) + facet_grid(.~dv, scales = "free_y", space="free")+  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +  geom_pointrange(aes(x = controls, colour=controls, y = X1, ymin = X1 - X2*interval2, ymax = X1 + X2*interval2), lwd = 1/2, position = position_dodge(width = 1/2), shape = 21) + theme_bw() + xlab("Control variables included (Yes/No)") + ylab("Unstandardized OLS regression coefficient")+ geom_linerange(aes(x = controls, colour=controls, ymin = X1 - X2*interval1,  ymax = X1 + X2*interval1), lwd = 1, position = position_dodge(width = 1/2))+ scale_color_manual(values=c('Grey','Black'), guide = FALSE)
ggsave(main, file="Main_results_noweights.pdf")

#Appendix Figure  A.4 plot of survey completes per day--------------------- 
plot_days<-ggplot(data, aes(x=day)) + geom_histogram(binwidth=0.5,colour = "black") + geom_vline(data=data, aes(xintercept = 9), colour="red") + theme_bw() + ylab("Completed surveys") + xlab("March")
ggsave(plot_days, file="Plot_days.pdf")


savehistory(file = "History.Rhistory")

